This is the analysis file corresponding to the manuscript “In search of the Goldilocks zone for hybrid speciation II: hard times for hybrid speciation”.
Default values to debug the add_dataframe() function.
if (FALSE){
k="output_sym_dom_rfunction_pi_0d5_rdot20.txt"
root="output_sym_adjabab_pi_0d5_nos"
}
Loading functions and packages
source("functions.R")
library(ggplot2)
## Warning: le package 'ggplot2' a été compilé avec la version R 4.2.2
print(paste(c("ggplot version ", as.character(packageVersion("ggplot2"))),collapse=""))
## [1] "ggplot version 3.4.0"
library(fields)
## Warning: le package 'fields' a été compilé avec la version R 4.2.2
## Le chargement a nécessité le package : spam
## Warning: le package 'spam' a été compilé avec la version R 4.2.2
## Spam version 2.9-1 (2022-08-07) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attachement du package : 'spam'
## Les objets suivants sont masqués depuis 'package:base':
##
## backsolve, forwardsolve
## Le chargement a nécessité le package : viridis
## Warning: le package 'viridis' a été compilé avec la version R 4.2.2
## Le chargement a nécessité le package : viridisLite
## Warning: le package 'viridisLite' a été compilé avec la version R 4.2.2
##
## Try help(fields) to get started.
print(paste(c("fields version ", as.character(packageVersion("fields"))),collapse=""))
## [1] "fields version 14.1"
library(viridis)
print(paste(c("viridis version ", as.character(packageVersion("viridis"))),collapse=""))
## [1] "viridis version 0.6.2"
library(RColorBrewer)
print(paste(c("RColorBrewer version ", as.character(packageVersion("RColorBrewer"))),collapse=""))
## [1] "RColorBrewer version 1.1.3"
library(cowplot)
## Warning: le package 'cowplot' a été compilé avec la version R 4.2.2
print(paste(c("cowplot version ", as.character(packageVersion("cowplot"))),collapse=""))
## [1] "cowplot version 1.1.1"
library(plyr)
## Warning: le package 'plyr' a été compilé avec la version R 4.2.2
print(paste(c("plyr version ", as.character(packageVersion("plyr"))),collapse=""))
## [1] "plyr version 1.8.8"
Load default dataset (1,000 replicates per combination of parameter)- Per default load the resulting dataset directly.
if (FALSE){
whole_dataset=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
colnames(whole_dataset)=names_col
known_folder=c()
} else{
whole_dataset=read.csv("whole_dataset.csv")
#whole_dataset$arch=as.factor(list_Arch[whole_dataset$arch])
known_folder=read.table("know_folder.txt")
}
Retrieve all possible folders starting by “output_” in the working directory.
list_folder=list.dirs(path = ".", full.names = TRUE, recursive = TRUE)
list_folder=gsub("./","",list_folder)
list_folder=list_folder[grep("^output", list_folder)]
Check the current main dataframe and update it if necessary
update=FALSE
for (i in list_folder){
if (!(i %in% unlist(known_folder))){
update=TRUE
#check_non_finish(i)
temp_dataset=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
temp_dataset2=add_dataframe(i,temp_dataset,truncature = 1000)
temp_dataset2=temp_dataset2[-c(1,2),]
colnames(temp_dataset2)=names_col
if(!check_arch(temp_dataset2)){
print(i)
break
}
whole_dataset=rbind(whole_dataset,temp_dataset2)
known_folder=c(known_folder,i)
}
if(length(known_folder)%%40==0){print(length(known_folder)/length(list_folder))}
}
whole_dataset=whole_dataset[!is.na(whole_dataset$n_rep),]
if (update){
write.table(known_folder,"know_folder.txt")
write.csv(whole_dataset,"whole_dataset.csv",row.names = FALSE)
}
Reload the data for better consistency
whole_dataset=read.csv("whole_dataset.csv")
whole_dataset$arch=as.factor(whole_dataset$arch)
whole_dataset=whole_dataset[-c(1,2),]
Load detailed dataset (10,000 replicates per combination of parameter) Per default load the resulting dataset directly.
if (FALSE){
whole_dataset_extra=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
colnames(whole_dataset_extra)=names_col
known_folder_extra=c()
dismiss_folder_extra=c()
} else{
whole_dataset_extra=read.csv("whole_dataset_extra.csv")
#whole_dataset_extra$arch=as.factor(list_Arch[whole_dataset_extra$arch])
known_folder_extra=read.table("know_folder_extra.txt")
}
Check the current main dataframe and update it if necessary
update=FALSE
temp_dataset=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
for (i in list_folder){
if (!(i %in% unlist(known_folder_extra))){
update=TRUE
temp_dataset2=add_dataframe(i,temp_dataset,truncature = 10000)
temp_dataset2=temp_dataset2[-c(1,2),]
colnames(temp_dataset2)=names_col
if(!check_arch(temp_dataset2)){
print(i)
break
}
whole_dataset_extra=rbind(whole_dataset_extra,temp_dataset2)
known_folder_extra=c(known_folder_extra,i)
}
if(length(known_folder_extra)%%40==0){print(length(known_folder_extra)/length(list_folder))}
}
whole_dataset_extra=whole_dataset_extra[!is.na(whole_dataset_extra$n_rep),]
if (update){
write.table(known_folder_extra,"know_folder_extra.txt")
write.csv(whole_dataset_extra,"whole_dataset_extra.csv",row.names = FALSE)
}
Reload the data for better consistency
whole_dataset_extra=read.csv("whole_dataset_extra.csv")
whole_dataset_extra$arch=as.factor(whole_dataset_extra$arch)
whole_dataset_extra=whole_dataset_extra[-c(1,2),]
Relationship between survival, epistasis and initial population size. (Since all architectures are equivalent, we always used the Adjacent ABAB architecture for the case of unlinked loci.)
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$r12==0.5&whole_dataset$arch %in% c("adj_ABAB","sin_DMI")&whole_dataset$K==1000000&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000)&whole_dataset$s1<0,]
#pdf("Figure/P_alive_vs_ep_unlinked.pdf",width=8, height=6)
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
temp_name2=c(`0`="One DMI",`1`="Two DMIs")
names(temp_name2)=c("0","1")
extra<-data.frame(label=c("A","B","C","D"),h12=c(0,1,0,1),arch=c("sin_DMI","sin_DMI","adj_ABAB","adj_ABAB"))
# adding extra column
temp_data$Cond=sapply(c(1:length(temp_data[,1])), function(x) paste(temp_data$h12[x], temp_data$arch[x], sep="_"))
temp_data$extra=plyr::mapvalues(temp_data$Cond, c("0_sin_DMI", "1_sin_DMI" , "0_adj_ABAB","1_adj_ABAB"), c("A","B","C","D"))
#str(temp_data)
ggplot(temp_data,aes(x=-e12,y=p_alive/n_rep,col=as.factor(N)))+geom_point()+scale_x_log10()+ xlab("Strength of epistasis")+ylab("Prob(alive)")+geom_line(size=0.5)+labs(col = "Initial pop. size")+geom_text(x=-4,y=0.9,aes(label=extra), colour="black")+facet_grid(as.numeric(temp_data$arch!="sin_DMI")~h12,labeller =labeller(.rows=temp_name2,.cols=temp_name))+theme(legend.position="bottom")+ylim(0,1)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
#dev.off()
Eeffect of epistasis on the resolution time of the genetic conflict of a single DMI and two DMIs.
#pdf("Figure/Res_time_vs_ep_unlinked.pdf",width=8, height=6)
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$r12==0.5&whole_dataset$arch %in% c("sin_DMI")&whole_dataset$K==1000000&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000)&whole_dataset$s1<0,]
ggplot(temp_data,aes(x=-e12,y=t_alive_dmi1,col=as.factor(N)))+geom_point()+scale_x_log10()+ xlab("Strength of epistasis")+ylab("Resolution time")+labs(col = "Initial pop. size")+facet_grid(as.numeric(temp_data$arch!="sin_DMI")~h12,labeller =labeller(.rows=temp_name2,.cols=temp_name))+theme(legend.position="bottom")+geom_line()+scale_y_log10()
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$r12==0.5&whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$K==1000000&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000)&whole_dataset$s1<0,]
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
temp_name2=c(`0`="One DMI",`1`="Two DMIs")
names(temp_name2)=c("0","1")
ggplot(temp_data,aes(x=-e12,y=t_alive_both,col=as.factor(N)))+geom_point()+scale_x_log10()+ xlab("Strength of epistasis")+ylab("Resolution time")+labs(col = "Initial pop. size")+facet_grid(as.numeric(temp_data$arch!="sin_DMI")~h12,labeller =labeller(.rows=temp_name2,.cols=temp_name))+theme(legend.position="bottom")+geom_line()+scale_y_log10()
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
#geom_smooth(se=F)+ xlim(-1,0)+
#dev.off()
Relationship between survival, epistasis and initial population size, as above, but with population size on the X axis.
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$K==1000000&whole_dataset$s1==-0.001&whole_dataset$r12==0.5&whole_dataset$arch %in% c("adj_ABAB","sin_DMI")&whole_dataset$e12 %in% c(0,-0.001,-0.01,-0.1,-0.5,-0.99),]
#pdf("Figure/P_alive_vs_Nini_unlinked.pdf",width=8, height=6)
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(-e12)))+geom_point(size=1)+geom_line()+scale_x_log10()+ xlab("Initial pop. size")+ylab("Prob(alive)")+labs(col = "Strength of epistasis")+facet_grid(as.numeric(temp_data$arch!="sin_DMI")~h12,labeller =labeller(.rows=temp_name2,.cols=temp_name))+theme(legend.position="bottom")+ylim(0,1)
#dev.off()
Relationship between minimum population size and epistasis
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$r12==0.5&whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000),]
#pdf("Figure/Min_pop_size_vs_ep_unlinked.pdf",width=8, height=6)
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
ggplot(temp_data,aes(x=-e12,y=N_min_alive,col=as.factor(N)))+geom_point()+geom_line()+ xlab("Strength of epistasis")+ylab("Minimum pop. size")+labs(col = "Initial pop. size")+xlim(0,1)+facet_wrap(~h12,labeller =as_labeller( temp_name))+theme(legend.position="bottom")+coord_cartesian(ylim = c(0, 500))+scale_x_log10()
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
#dev.off()
Relationship between resolution time and epistasis
#pdf("Figure/Time_resolution_vs_ep_unlinked.pdf",width=8, height=6)
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
ggplot(temp_data,aes(x=-e12,y=t_both,col=as.factor(N)))+geom_point() +geom_line()+ xlab("Strength of epistasis")+ylab("Time resolution 2 DMIs")+labs(col = "Initial pop. size")+xlim(0,1)+facet_wrap(~h12,labeller =as_labeller( temp_name))+theme(legend.position="bottom")+coord_cartesian(ylim = c(0, 100))
#dev.off()
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$r12==0.5&whole_dataset$h12==1 &whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$N %in% c(100,500,1000,2500,5000,10000),] # codominant case
temp_data2=whole_dataset[default_regime(whole_dataset)&whole_dataset$r12==0.5&whole_dataset$h12==0 &whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$N %in% c(100,500,1000,2500,5000,10000),]#recessive case
ggplot(temp_data,aes(x=N_min_alive,y=t_both,col=p_alive/n_rep))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution 2 codominant DMIs")+labs(col = "P(alive)")+theme(legend.position="bottom")+scale_x_log10()+scale_y_log10()
ggplot(temp_data,aes(x=N_min_alive,y=t_both,col=as.factor(N)))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution 2 codominant DMIs")+labs(col = "Initial pop. size")+theme(legend.position="bottom")+scale_x_log10()+scale_y_log10()
ggplot(temp_data2,aes(x=N_min_alive,y=t_both,col=p_alive/n_rep))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution 2 recessive DMIs")+labs(col = "P(alive)")+theme(legend.position="bottom")+scale_x_log10()+scale_y_log10()
ggplot(temp_data2,aes(x=N_min_alive,y=t_both,col=as.factor(N)))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution 2 recessive DMIs")+labs(col = "Initial pop. size")+theme(legend.position="bottom")+scale_x_log10()+scale_y_log10()
To better understand the outcome observed above, and why intermediate epistasis values presented the highest risk of extinction, we focused on the deterministic scenario. ### For an initial pop size of N=5000
data0=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep0.txt")
data1=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.1.txt")
data2=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.2.txt")
data3=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.3.txt")
data4=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.4.txt")
data5=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.5.txt")
data6=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.6.txt")
data7=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.7.txt")
data8=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.8.txt")
data9=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.9.txt")
data99=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.99.txt")
Population size evoltuion for different epistasis strengths
#pdf("Figure/Determ_pop_size_N5000_unlinked.pdf",height=4.5,width=6)
plot(data5$V1,type='l',ylim=c(4,100),xlim=c(0,200),col="purple",xlab="Generations",ylab="Pop. size",lwd=2,log="y")
points(seq(1000),data6$V1,type='l',col="blue",lwd=2)
points(seq(1000),data7$V1,type='l',col="cyan",lwd=2)
points(seq(1000),data8$V1,type='l',col="orange",lwd=2)
points(seq(1000),data9$V1,type='l',col="red",lwd=2)
legend("topright",legend=c(expression(paste(epsilon,"=-0.5")),expression(paste(epsilon,"=-0.6")),expression(paste(epsilon,"=-0.7")),expression(paste(epsilon,"=-0.8")),expression(paste(epsilon,"=-0.9"))),col=c("purple","blue","cyan","orange","red"),pch=16)
#dev.off()
We want to calculate here the initial mean growth rate of the population to compare it to en evolutionary rescue scenario (which selection coefficient in a single locus model will generate the same initial and final growth rate).
#pdf("Figure/rbar-vs-ep.pdf",width=8,height=6)
plot(-c(0,0.1,0.2,0.3,0.4,.5,.6,.7,.8,.9,.99),c(data0$V1[2]/data0$V1[1],data1$V1[2]/data1$V1[1],data2$V1[2]/data2$V1[1],data3$V1[2]/data3$V1[1],data4$V1[2]/data4$V1[1],data5$V1[2]/data5$V1[1],data6$V1[2]/data6$V1[1],data7$V1[2]/data7$V1[1],data8$V1[2]/data8$V1[1],data9$V1[2]/data9$V1[1],data99$V1[2]/data99$V1[1]),xlab=expression(paste("Epistasis, ",epsilon)),ylab="Average r",pch=16,col="blue",type="b",cex.axis=1.2,cex.lab=1.2,ylim=c(0.4,1.01))
legend("topleft",legend = c("F1","F2"),col=c("blue","cyan"),pch=16,cex=1.2)
points(-c(0,0.1,0.2,0.3,0.4,.5,.6,.7,.8,.9,.99),c(data0$V1[2]/data0$V1[1],data1$V1[3]/data1$V1[2],data2$V1[3]/data2$V1[2],data3$V1[3]/data3$V1[2],data4$V1[3]/data4$V1[2],data5$V1[3]/data5$V1[2],data6$V1[3]/data6$V1[2],data7$V1[3]/data7$V1[2],data8$V1[3]/data8$V1[2],data9$V1[3]/data9$V1[2],data99$V1[3]/data99$V1[2]),col="cyan",pch=16,type="b")
#dev.off()
c(data0$V1[2]/data0$V1[1],data1$V1[2]/data1$V1[1],data2$V1[2]/data2$V1[1],data3$V1[2]/data3$V1[1],data4$V1[2]/data4$V1[1],data5$V1[2]/data5$V1[1],data6$V1[2]/data6$V1[1],data7$V1[2]/data7$V1[1],data8$V1[2]/data8$V1[1],data9$V1[2]/data9$V1[1],data99$V1[2]/data99$V1[1])
## [1] 1.0059661 0.9103993 0.8248922 0.7494448 0.6840570 0.6287288 0.5834604
## [8] 0.5482516 0.5231024 0.5080129 0.5030334
sqrt(c(data1$V1[2]/data1$V1[1],data2$V1[2]/data2$V1[1],data3$V1[2]/data3$V1[1],data4$V1[2]/data4$V1[1],data5$V1[2]/data5$V1[1],data6$V1[2]/data6$V1[1],data7$V1[2]/data7$V1[1],data8$V1[2]/data8$V1[1],data9$V1[2]/data9$V1[1],data99$V1[2]/data99$V1[1])/1.005966)-1
## [1] -0.04868509 -0.09446145 -0.13686613 -0.17537882 -0.20943053 -0.23842262
## [7] -0.26175878 -0.27888967 -0.28936640 -0.29285779
#pdf("Figure/equivalent_s.pdf",width=8,height=6)
plot(-c(0,0.1,0.2,0.3,0.4,.5,.6,.7,.8,.9,.99),sqrt(c(1,data1$V1[2]/data1$V1[1],data2$V1[2]/data2$V1[1],data3$V1[2]/data3$V1[1],data4$V1[2]/data4$V1[1],data5$V1[2]/data5$V1[1],data6$V1[2]/data6$V1[1],data7$V1[2]/data7$V1[1],data8$V1[2]/data8$V1[1],data9$V1[2]/data9$V1[1],data99$V1[2]/data99$V1[1])/1.005966)-1,pch=16,xlab="Epistasis",ylab="Equivalent s")
#dev.off()
Same as above but with a reduction by a factor 10 of the initial population size.
data3=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.3.txt")
data4=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.4.txt")
data5=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.5.txt")
data6=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.6.txt")
data7=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.7.txt")
data8=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.8.txt")
data9=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.9.txt")
data99=read.table("output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.99.txt")
#pdf("Figure/Determ_pop_size_N500_unlinked.pdf",height=4.5,width=6)
plot(data5$V1,type='l',ylim=c(3,500),xlim=c(0,750),col="purple",xlab="Generations",ylab="Pop. size",lwd=2,log="y")
points(seq(1000),data6$V1,type='l',col="blue",lwd=2)
points(seq(1000),data7$V1,type='l',col="cyan",lwd=2)
points(seq(1000),data4$V1,type='l',col="darkorchid4",lwd=2)
points(seq(1000),data3$V1,type='l',col="black",lwd=2)
legend("bottomright",legend=c(expression(paste(epsilon,"=-0.3")),expression(paste(epsilon,"=-0.4")),expression(paste(epsilon,"=-0.5")),expression(paste(epsilon,"=-0.6")),expression(paste(epsilon,"=-0.7"))),col=c("black","darkorchid4","purple","blue","cyan","orange"),pch=16)
#dev.off()
#Analysis of the linked loci: Codominant scenario
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$h12==1&whole_dataset$arch!="sin_DMI",]
Legned for the 6 architectures
#pdf("Figure/legend_arch.pdf",width=7,height=1)
plot_grid(get_legend(ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point(size=4)+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20),legend.position = "bottom")+ylim(0,1)))
#dev.off()
Relationship between survival probability and minimum population size for all 6 architectures
#pdf("Figure/Min_pop_size_vs_alive_cod.pdf",width=8,height=6)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=r12))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Recombination rate")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
ggplot(temp_data,aes(x=t_both,y=p_alive/n_rep,col=r12))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the 2 DMIs")+ylab("Prob(alive)")+labs(col = "Recombination rate")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=log10(t_both)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "(log10) Time of resolution of the 2 DMIs")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
ggplot(temp_data,aes(x=N,y=N_min_alive,col=p_alive/n_rep))+geom_point()+geom_jitter(width=0.2)+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Starting population size")+ylab("Minimum population size")+labs(col = "P(alive)")+scale_x_log10()+scale_y_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom")
ggplot(temp_data,aes(x=t_both,y=N_min_alive,col=p_alive/n_rep))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the 2 DMIs")+ylab("Minimum population size")+labs(col = "P(alive)")+scale_x_log10()+scale_y_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
ggplot(temp_data,aes(x=t_both,y=p_alive/n_rep,col=log10(N_min_alive)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the 2 DMIs")+ylab("Prob(alive)")+labs(col = "Minimum population size")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
#dev.off()
Same as above but with guidelines
#pdf("Figure/Min_pop_size_vs_alive_cod_guide.pdf",width=8,height=6)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=log10(t_both)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("Prob(survival)")+labs(col = "(log10) Time of resolution of the 2 DMIs")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))+ geom_segment(aes(x = 10, y = 0, xend = 1000, yend = 1),size=1,color="black")
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=(r12)))+geom_point(size=0.5)+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("Prob(survival)")+labs(col = "Recombination rate")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))+ geom_segment(aes(x = 10, y = 0, xend = 600, yend = 1),size=1,color="black")+ theme(axis.text=element_text(size=10),axis.title =element_text(size=12) )+theme(strip.text.x = element_text(size = 12))+guides(color=FALSE)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
#dev.off()
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$h12==1&whole_dataset$arch!="sin_DMI"&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000),]
#pdf("Figure/p_alive_vs_Nmin.pdf",width=8, height=6)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=as.factor(N)))+geom_point()+ xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Initial Population size")+scale_x_log10()+ theme(legend.position="bottom")+scale_colour_brewer(palette="Set1")+ylim(0,1)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+ xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Initial Population size")+scale_x_log10()+ theme(legend.position="bottom")+scale_colour_brewer(palette="Set1")+ylim(0,1)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=as.factor(N)))+geom_point()+ xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Initial Population size")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom")+scale_colour_brewer(palette="Set1")+ylim(0,1)
#dev.off()
Relationship between hybrid speciation and recombination for all 6 architectures and different population sizes
#pdf("all_6arch_hard_sel.pdf",width=8,height=6)
#adj ABAB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Adjacent ABAB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#adj ABBA
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Adjacent ABBA DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#cro AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Crosed ABBA DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#cro AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Crossed AABB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#nes ABAB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Nested ABAB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#nes AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Nested AABB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#dev.off()
Relationship between hybrid speciation conditioned on survival and recombination for all 6 architectures and different population sizes
#pdf("all_6arch_hard_sel_cond.pdf",width=8,height=6)
#adj ABAB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Adjacent ABAB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#adj ABBA
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Adjacent ABBA DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#cro AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Crosed ABBA DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#cro AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Crossed AABB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#nes ABAB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Nested ABAB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#nes AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Nested AABB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#dev.off()
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==5000&whole_dataset$arch!="sin_DMI",]
Relationship between survival probability, hybrid speciation and recombination for all architectures
#pdf("Figure/cod_case_ep0d2_N5000.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylim(0.4,1)+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.25) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid sp.|survival)")+xlab("Recombination")+geom_line()#+ guides(col=FALSE)
## Warning: Removed 1 rows containing missing values (`geom_point()`).
#dev.off()
relationship between minimum population size, hybrid speciation probbaility and recombination
#pdf("Figure/Nminvsrec_and_hyybvsrec_cod0d2.pdf",widt=8,height=6)
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Min. pop. size")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
#dev.off()
Further exploration of this scenario, including time of resolution of the DMIs
ggplot(temp_data,aes(y=t_alive_both,x=N_min_alive,col=p_alive/n_rep))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution")+labs(col = "P(alive)")+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom")+scale_colour_gradient(low = "orange", high = "purple")
ggplot(temp_data,aes(x=N_alive_both,y=t_alive_both,col=p_alive/n_rep))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Population size when the DMIs are resolved")+ylab("Time of resolution of the two DMIs")+labs(col = "Prob(surviving)")
ggplot(temp_data,aes(x=N_alive_both,y=t_alive_both,col=((p_hybA1+p_hybB1)/p_alive)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Population size when the DMIs are resolved")+ylab("Time of resolution of the two DMIs")+labs(col = "Prob(hyb_spec|alive)")
ggplot(temp_data,aes(x=N_alive_both,y=t_alive_both,col=r12))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Population size when the DMIs are resolved")+ylab("Time of resolution of the two DMIs")+labs(col = "Recombination rate")
ggplot(temp_data,aes(x=N_alive_both,y=p_alive/n_rep,col=r12))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Population size when the DMIs are resolved")+ylab("Prob(alive)")+labs(col = "Recombination rate")
ggplot(temp_data,aes(x=t_alive_both,y=p_alive/n_rep,col=log10(r12),pch=as.factor(arch)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the two DMIs")+ylab("Prob(alive)")+labs(col = "log(Recombination rate)")
ggplot(temp_data,aes(x=p_alive/n_rep,y=(p_hybA1+p_hybB1)/n_rep,col=r12,pch=as.factor(arch)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+ylab("Prob(Hybrid speciation)")+xlab("Prob(alive)")+labs(col = "log(Recombination rate)")
ggplot(temp_data,aes(x=p_alive/n_rep,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=r12,pch=as.factor(arch)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+ylab("Prob(Hybrid speciation|alive)")+xlab("Prob(alive)")+labs(col = "log(Recombination rate)")
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=r12))+geom_point()+ xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Recombination")+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom")+scale_colour_gradient(low = "orange", high = "purple")
Here we use a more precise estimation of the different outcomes
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==5000&whole_dataset_extra$arch!="sin_DMI",]
Relationship between hybrid pseciation, survival probability ansd recombination, with higher resolution
#pdf("Figure/cod_case_ep0d2_N5000_10k.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylim(0.4,1)+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.25) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid sp.|survival)")+xlab("Recombination")+geom_line()#+ guides(col=FALSE)
#dev.off()
Relationship between minimum population size and hybrid speciation with higher resolution
#pdf("Figure/Nminvsrec_and_hyybvsrec_cod0d2_10k.pdf",widt=8,height=6)
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10(breaks=c(5*10^-4, 10^-2,0.5))+ylab("Min. pop. size")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(axis.text.x = element_text(angle=45),plot.margin=margin(5.5,15,0,0,"pt"))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label)) + theme(axis.text.x = element_text(angle=45),plot.margin=margin(5.5,15,0,0,"pt"))
#dev.off()
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==500,]
#pdf("Figure/cod_case_ep0d2_N500.pdf",widt=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylim(0,1)+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.1)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.15) + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()
#dev.off()
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$K==1000000&whole_dataset_extra$arch!="sin_DMI"&whole_dataset_extra$s1==-0.001&whole_dataset_extra$s2==-0.001&whole_dataset_extra$s3==-0.001&whole_dataset_extra$s4==-0.001,]
#pdf("Figure/cod_case_ep0d2_N100.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()
## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
#dev.off()
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==50000,]
Relation between survival probability and hybrid speciation in large population
#pdf("Figure/cod_case_ep0d2_N50000.pdf",widt=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylim(0,1)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10() + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")
#dev.off()
For the Adjacent ABAB architecutre
temp_data=whole_dataset[whole_dataset$arch=="adj_ABAB"&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2)&whole_dataset$h12==1&whole_dataset$N %in% c(500,5000,50000),]
#pdf("Figure/combined_surv_hyb_sp.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(N)))+geom_point()+geom_point(aes(x=r12,y=p_alive/n_rep,col=as.factor(N)),shape=0)+xlab("Recombination")+ylab("P(hybrid sp.) or P(alive)")+ scale_colour_manual( values=seq(8), name="Initial pop size") +theme(text = element_text(size=20))+ylim(0,1)+scale_x_log10()+theme(legend.position="bottom")
#ggplot(temp_data,aes(x=p_alive/n_rep,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch),labeller = as_labeller(arch_label)))+geom_point()+xlab("P(alive)")+ylab("P(hybrid sp.)")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylim(0,1)+ guides(col=FALSE)
#dev.off()
###Investigating the effect of direct selection (or its lack) Same as above, but we are removing the effect of direct selection on all four loci
temp_data=whole_dataset_extra[whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$K==1000000&whole_dataset_extra$arch!="sin_DMI"&whole_dataset_extra$s1==0&whole_dataset_extra$s2==0&whole_dataset_extra$s3==0&whole_dataset_extra$s4==0,]
#pdf("Figure/cod_case_ep0d2_N100_10k_nosel.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()
## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
#dev.off()
temp_data=whole_dataset_extra[whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$arch!="sin_DMI"&whole_dataset_extra$K==1000000,]
Comparaison between the case with and without direct selection
#pdf("Figure/comaprison_no_sel.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(s1)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination") +theme(text = element_text(size=20))+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Selection")#+ guides(col=FALSE)#+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture")
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(s1)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination") +theme(text = element_text(size=20))+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Selection")
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(s1)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=20))+ylab("P(hybrid speciation&survival)")+xlab("Recombination")+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Selection")
## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(s1)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Selection")
## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
# dev.off()
temp_data=whole_dataset[whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==100&whole_dataset$arch!="sin_DMI"&whole_dataset$s1==-0.001&whole_dataset$s2==-0.001&whole_dataset$s3==-0.001&whole_dataset$s4==-0.001,]
Same as above, but we are removing the varying the carrying capacity of the environment
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(K)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination") +theme(text = element_text(size=20))+facet_wrap(~arch)#+ guides(col=FALSE)#+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture")
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(K)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination") +theme(text = element_text(size=20))+facet_wrap(~arch)#+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture")
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(K)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=20))+ylab("P(hybrid speciation&survival)")+xlab("Recombination")+facet_wrap(~arch)+geom_smooth(se=F)+ guides(col=FALSE)
## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(K)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+facet_wrap(~arch)+geom_smooth(se=F)+ guides(col=FALSE)#+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture")
## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))
# ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))
# ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_point()
# ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(arch)))+geom_point(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation&survival)")+xlab("Recombination")
# ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(arch)))+geom_point(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation&survival)")+xlab("Recombination")+geom_smooth(se=F)
Same, but with higher resolution
temp_data=whole_dataset_extra[whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$arch!="sin_DMI"&whole_dataset_extra$s1==-0.001&whole_dataset_extra$s2==-0.001&whole_dataset_extra$s3==-0.001&whole_dataset_extra$s4==-0.001,]
#pdf("Figure/as_a_fucntion_of_K.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(K)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination") +theme(text = element_text(size=20))+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Max pop. size")#+ guides(col=FALSE)#+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture")
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(K)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination") +theme(text = element_text(size=20))+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Max pop. size")
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(K)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=20))+ylab("P(hybrid speciation&survival)")+xlab("Recombination")+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Max pop. size")
## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(K)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Max pop. size")
## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
#dev.off()
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==100000&whole_dataset$r12==0.1&whole_dataset$s1==-0.001&whole_dataset$s2==-0.001&whole_dataset$s3==-0.001&whole_dataset$s4==-0.001,]
Effect of the mean offspring number on survival and hybrid specation in alrge populations
ggplot(temp_data,aes(x=mean_off,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("Mean off")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+xlim(1,1.03)
ggplot(temp_data,aes(x=mean_off,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(survival)")+xlab("Mean off")+xlim(1,1.015)+geom_line()
## Warning: Removed 89 rows containing missing values (`geom_point()`).
## Warning: Removed 89 rows containing missing values (`geom_line()`).
ggplot(temp_data,aes(x=mean_off,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Mean off")+xlim(1,1.03)
## Warning: Removed 5 rows containing missing values (`geom_point()`).
Effect of the mean offspring number on survival and hybrid specation in large populations
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==1000&whole_dataset$r12==0.1&whole_dataset$s1==-0.001&whole_dataset$s2==-0.001&whole_dataset$s2==-0.001&whole_dataset$s4==-0.001,]
ggplot(temp_data,aes(x=mean_off,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("Mean off")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))
ggplot(temp_data,aes(x=mean_off,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(survival)")+xlab("Mean off")
ggplot(temp_data,aes(x=mean_off,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Mean off")
## Warning: Removed 5 rows containing missing values (`geom_point()`).
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$K==1000000&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$r12==0.1&whole_dataset$s1==-0.001&whole_dataset$s2==-0.001&whole_dataset$s2==-0.001&whole_dataset$s4==-0.001,]
Relationship between survival and hybrid speciation and the initial population size
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(survival)")+xlab("N_ini")+scale_x_log10()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("N_ini")+scale_x_log10()
Zoom on the lower values for survival
#pdf("Figure/surv_vs_Nin_cod_ep0d2_r0d1.pdf",width=8,height=6)
ggplot(temp_data[temp_data$N<=1000,],aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(survival)")+xlab("N_ini")+scale_x_log10()+geom_line()
#dev.off()
With higher resolution
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$s1==-0.001&whole_dataset_extra$s2==-0.001&whole_dataset_extra$s3==-0.001&whole_dataset_extra$s4==-0.001&whole_dataset_extra$K==1000000,]
#pdf("Figure/hybrid_spec_Nini.pdf",width=8, height=6)
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_line() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("N_ini")+scale_x_log10()+geom_point()
#dev.off()
Decomposition of the outcome as a function of the initial population size
ggplot(temp_data,aes(x=N,y=(p_alive)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(alive)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()+ylim(0.05,.2)+scale_x_log10(limits=c(10,750))
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 74 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 74 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation)")+xlab("N_ini")+ scale_colour_manual(values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation)")+xlab("N_ini")+ scale_colour_manual(values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_point()+ylim(0,0.05)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 35 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 35 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation & surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()+ylim(0,0.05)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 10 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 10 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A1B2 & surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A2B1 & surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_alive_both,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A2B1 & surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_alive_both,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Time resolution of DMIs")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()+scale_y_log10()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_diff,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Difference in resolution of DMIs")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_diff_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Difference in resolution of DMIs|alive")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_diff_alive_hybA1,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Difference in resolution of DMIs|alive for hyb A1")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=t_diff_alive_hybB1,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Difference in resolution of DMIs|alive for hyb B1")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(Alive)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid sp.|surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_line()+ylab("P(hybrid sp.|surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation|surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A1B2|surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()+ylim(0,0.05)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 1 rows containing missing values (`geom_smooth()`).
ggplot(temp_data,aes(x=N,y=(p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A2B1|surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()+ylim(0,0.05)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 20 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 20 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+ylab("P(Alive)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()+geom_line()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+ylab("P(hybrid speciation|surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()+geom_smooth(se=F)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A1B2|surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A2B1|surviving)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_point()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Stacks plot representing the different outcomes for the different
architectures
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="adj_ABBA"&!is.na(whole_dataset_extra$arch)&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_adjABBA.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Adjacent ABBA")
#dev.off()
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="adj_ABAB"&!is.na(whole_dataset_extra$arch)&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_adjABAB.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Adjacent ABAB")
#dev.off()
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="cro_AABB"&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_croAABB.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Crossed AABB")
#dev.off()
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="cro_ABBA"&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_croABBA.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Crossed ABBA")
#dev.off()
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="nes_AABB"&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_nesAABB.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Nested AABB")
#dev.off()
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="nes_ABAB",]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_nesABAB.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Nested ABAB")
#dev.off()
Decomposition with higher resolution
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01,]
ggplot(temp_data,aes(y=(p_alive_hybA1+p_alive_hybB1)/n_rep,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_y_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(Hybrid speciation&survival)")+xlab("Population size at the time of resolution of the genetic conflicts")
#ggplot(temp_data,aes(y=NL_alive_A1,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_y_log10()
#ggplot(temp_data,aes(y=NL_alive_A2,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_y_log10()
#ggplot(temp_data,aes(y=NL_alive_B1,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_y_log10()
#ggplot(temp_data,aes(y=NL_alive_B2,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_y_log10()
#temp_soft=read.table("~/IGC/dmi_and_pattern/temp_soft_sel.csv",sep=",",header=TRUE)
#ggplot(temp_soft,aes(x=N,y=(p_hybA1 +p_hybB1)/1000,col=arch ))+geom_point()+theme(text = element_text(size=20))+scale_x_log10()+ylim(0,0.2)+ylab("Hybdrid spe.")+geom_line()
Relation between survival and initial population size for unlinked loci
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$r12==0.5&whole_dataset$mean_off==1.01,]
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(survival)")+xlab("N_ini")+scale_x_log10()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("N_ini")+scale_x_log10()
Relation between survival and initial population size for closely linked loci
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$r12==0.01&whole_dataset$mean_off==1.01,]
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(survival)")+xlab("N_ini")+scale_x_log10()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("N_ini")+scale_x_log10()
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.99|whole_dataset$e13==-0.99|whole_dataset$e14==-0.99)&whole_dataset$h12==1&whole_dataset$N==5000&whole_dataset$arch!="sin_DMI",]
Relation between survival, hybrid speciation and recombination for (almost) lethal epistasis
#pdf("Figure/cod_case_ep0d99_N5000.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=(p_alive)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("P(survival)")+xlab("Recombination rate")+scale_colour_manual( values=seq(7), labels =list_Arch_legend ,name="Linkage\n architecture") +theme(text = element_text(size=20)) #+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.25) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+ guides(col=FALSE)
#dev.off()
Relation between survival and minimum population size for (almost) lethal epistasis
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Min. pop. size")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
####N=50000
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.99|whole_dataset$e13==-0.99|whole_dataset$e14==-0.99)&whole_dataset$h12==1&whole_dataset$N==50000,]
Relation between survival, hybrid speciation and recombination for (almost) lethal epistasis in large populations
ggplot(temp_data,aes(x=r12,y=(p_alive)/n_rep,col=as.factor(arch)))+geom_point()+geom_smooth()+scale_x_log10()+ylim(0,1)+ylab("Survival probability")+xlab("Recombination rate")+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.99|whole_dataset_extra$e13==-0.99|whole_dataset_extra$e14==-0.99)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$arch!="sin_DMI",]
Relation between survival, hybrid speciation and recombination for (almost) lethal epistasis in small populations with higher resolution
#pdf("Figure/cod_case_ep0d99_N100.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))
#dev.off()
#Analysis of the linked loci: Recessive scenario
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&whole_dataset$h12==0&whole_dataset$arch!="sin_DMI",]
Relationship between survival probability and minimum population size
#ggplot(temp_data,aes(x=t_alive_both,y=p_alive/n_rep,col=log10(r12),pch=as.factor(arch),labeller = as_labeller(arch_label)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the two DMIs")+ylab("Prob(alive)")+labs(col = "log(Recombination rate)")+scale_x_log10()
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=r12))+geom_point(size=0.5)+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("P(survival)")+labs(col = "Recombination rate")+scale_x_log10(breaks=c(10,1000,100000))+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.23, "in"))+geom_segment(aes(x = 10, y = 0, xend = 600, yend = 1),size=1,color="black")
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==0&whole_dataset$r12==0.1&whole_dataset$mean_off==1.01,]
Relationship between survival, hybrid speciation and initial population size for all 6 architectures
#pdf("Figure/rec_case_all_arch_vs_Nini.pdf",width=8,height=6)
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(survival)")+xlab("N_ini")+scale_x_log10()+geom_line()
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("N_ini")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+scale_x_log10()+geom_line()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("N_ini")+scale_x_log10()+geom_line()
#dev.off()
Relationship between survival, hybrid speciation and recombination for all 6 architectures
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==0&whole_dataset$N==5000&whole_dataset$arch!="sin_DMI",]
#pdf("Figure/rec_case_ep0d2_N5000.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylim(0.4,1)+geom_line()# guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()# guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()# guides(col=FALSE)
#dev.off()
Same as above with higher resolution
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==0&whole_dataset_extra$N==5000&whole_dataset_extra$arch!="sin_DMI",]
#pdf("Figure/rec_case_ep0d2_N5000_10k.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylim(0.85,1)+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()# guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+ geom_line()#guides(col=FALSE)
#dev.off()
Relationship between recombination, minimum population size and hybrid speciation, with higher resolution
#pdf("Figure/Nminvsrec_and_hyybvsrec_rec0d2_10k.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10(breaks=c(5*10^-4, 10^-2,0.5))+ylab("Min. pop. size")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(axis.text.x = element_text(angle=45),plot.margin=margin(5.5,15,0,0,"pt"))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10(breaks=c(500,1000, 2000,4000))+ylim(0,0.65)+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label)) + theme(axis.text.x = element_text(angle=45),plot.margin=margin(5.5,15,0,0,"pt"))
#dev.off()
###N=5000
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.99|whole_dataset$e13==-0.99|whole_dataset$e14==-0.99)&whole_dataset$h12==0&whole_dataset$N==5000&whole_dataset$arch!="sin_DMI",]
Relationship between survival, hybrid speciation and recombination for all 6 architectures for lethal epistasis
#pdf("Figure/rec_case_ep0d99_N5000.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylim(0.85,1)+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ geom_line()#guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()# guides(col=FALSE)
#dev.off()
###N=500
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.99|whole_dataset$e13==-0.99|whole_dataset$e14==-0.99)&whole_dataset$h12==0&whole_dataset$N==500&whole_dataset$arch!="sin_DMI",]
Relationship between survival, hybrid speciation and recombination for all 6 architectures for lethal epistasis and small populatiobs
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ geom_line()#guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()# guides(col=FALSE)
#Evolutionary rescue
Here we simulated with the same program an evolutionary rescue scenario, and evaluated the link between minimum population size and survival probability.
if (FALSE){
fake_dataset=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
colnames(fake_dataset)=names_col
} else{
fake_dataset=read.csv("fake_dataset.csv")
fake_dataset$arch=as.factor(list_Arch[fake_dataset$arch])
}
known_folder_f=read.table("know_folder_f.txt")
list_folder_f=list.dirs(path = ".", full.names = TRUE, recursive = TRUE)
list_folder_f=gsub("./","",list_folder_f)
list_folder_f=list_folder_f[grep("^fake", list_folder_f)]
update=FALSE
#known_folder_f=c()
for (i in list_folder_f){
if (!(i %in% unlist(known_folder_f))){
update=TRUE
#check_non_finish(i)
fake_dataset=add_dataframe(i,fake_dataset)
known_folder_f=c(known_folder_f,i)
}
}
#whole_dataset=whole_dataset[!is.na(whole_dataset$n_rep),]
if (update){
write.table(list_folder_f,"know_folder_f.txt")
write.csv(fake_dataset,"fake_dataset.csv",row.names = FALSE)
}
fake_dataset=read.csv("fake_dataset.csv")
fake_dataset$arch=as.factor(list_Arch[fake_dataset$arch])
fake_dataset=fake_dataset[-c(1,2),]
fake_dataset$N_min_alive[is.na(fake_dataset$N_min_alive)]=0
##Analysis
Link between survival (or evolutionary rescue) and minimum population size
#pdf("Figure/prob_rescue.pdf",width=8,height=6)
ggplot(fake_dataset,aes(x=N_min_alive,y=p_alive/n_rep,col=as.factor(f1000/N)))+geom_point()+scale_x_log10(limits=c(1,10^5))+ geom_segment(aes(x = 10, y = 0, xend = 600, yend = 1),size=1,color="black")+ylab("Prob. rescue")+xlab("Minimum population size")+labs(col="Freq. of\n allele A")
## Warning: Transformation introduced infinite values in continuous x-axis
#dev.off()
Link between strength of selection and minimum population size
ggplot(fake_dataset,aes(x=s1,y=N_min_alive,col=as.factor(f1000/N)))+geom_point()+ylab("Min pop size")+xlab("Selection on the ancestral allele")+labs(col="Freq. of allele A")+geom_line()